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Abstract. In this work we present a numerical method for the Optimal 
Mass Transportation problem. Optimal Mass Transportation (OT) is 
an active research field in mathematics. It has recently led to significant 
theoretical results as well as applications in diverse areas. Numerical 
solution techniques for the OT problem remain underdeveloped. The 
solution is obtained by solving the second boundary value problem for 
the Monge- Ampere equation, a fully nonlinear elliptic partial differential 
equation (PDE). Instead of standard boundary conditions the problem 
has global state constraints. These are reformulated as a tractable local 
PDE. We give a proof of convergence of the numerical method, using 
the theory of viscosity solutions. Details of the implementation and a 
fast solution method are provided in the companion paper [BFOI2]. 
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1. Introduction 



The Optimal Transportation (OT) problem is a simply posed mathe- 
matical problem which dates back more than two centuries. It has led to 
significant results in probability, analysis, and Partial Differential Equa- 
tions (PDEs), among other areas. The core theory is by now well estab- 
lished: good presentations are available the textbook [Vil03] and the sur- 
vey [Eva99]. The subject continues to find new relevance, in particular to 
geometry and curvature (see [Vil09]), and to dissipative evolutionary equa- 
tions (see [JK098, OttOl] and the resulting literature). Application areas 
include image registration, mesh generation, reflector design, astrophysics, 
and meteorology (see [Vil03] for references). 

This article presents a numerical method for the Optimal Transportation 
(OT) problem with quadratic costs. It is formally equivalent to the second 
boundary value problem for the Monge-Ampere equation, a fully nonlin- 
ear elliptic partial differential equation (PDE). In this setting the equation 
lacks standard boundary conditions: instead there are state constraints on 
the gradient of the solution. Because these state constraints are not easily 
enforced, we reformulate them using the distance function to the target set, 
which results in a local but implicit Hamilton-Jacobi PDE at the boundary. 
Our main result is a proof of convergence of solutions of the finite difference 
approximation to the PDE, using the theory of viscosity solutions. 

Going from a theoretically convergent algorithm to the simplest and most 
efficient numerical implementation requires more explanation, a different 
emphasis from the mathematical aspects of the convergence proof. For 
this reason, we present the details of the implementation in a companion 
paper, [BF012]. Basic computational results are presented here as an il- 
lustration. Additional computational results, including some which are not 
covered by the theory (e.g. non-convex targets, Alexsandroff solutions), can 
be found in [BF012]. 

1.1. Optimal Transportation and the Monge-Ampere PDE. The 

OT problem is described as follows. Suppose we are given two probabil- 
ity densities, where 



(1) 



Py(u) is Lipschitz continuous with convex support Y 
px{x) is bounded with support X, 
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and X,Y C W 1 are bounded and open. Consider the set, M, of maps which 
rearrange the measure py into the measure px, 

(2) M = {T : X i-> Y, p Y {T) det(VT) = p x }- 
The OT problem, in the case of quadratic costs, is given by 

(3) inf - / \\x — T(x)\\ 2 px(x)dx. 

See Figure 5 for an illustration of the optimal map between ellipses. 

General conditions under which the OT problem (2, 3) is well-posed are 
established in [Brc91]. The unique minimizing map, T, at which the mini- 
mum is reached, is the gradient of a convex function 

T = Vu, u convex iIcK^l, 

which is therefore also unique up to a constant. Write D 2 u for the Hessian 
of the function u. Formally substituting T = Vu into (2) results in the 
Monge- Ampere PDE 

(MA) det(D 2 u(x)) = P ^ X ] , for 

along with the restriction 

(C) u is convex. 

The PDE (MA) lacks standard boundary conditions. However, it is con- 
strained by the fact that the gradient map takes X to Y, 

(BV2) Vu(X) = Y. 

The condition (BV2) is referred to as the second boundary value problem for 
the Monge- Ampere equation in the literature (see [Urb97]). We will also use 
the term OT boundary conditions. 



Existing approaches to solving the OT problem using PDE methods fail 
to address the condition (BV2) for general densities. Instead, they generally 
limit their applicability to rectangular density supports and use either pe- 
riodic or inhomogeneous Neumann boundary conditions, provided the mass 
of the target density does not vanish. 

Another way to around this problem is to allow densities to vanish, how- 
ever this introduces new problems: Allowing px to vanish causes the the 
PDE (MA) to become degenerate elliptic. Allowing py to be discontinuous 
causes instabilities, since the Lipschitz constant of py appears in gradient 
descent type of methods. Smoothing and adding mass to bound either den- 
sity away from zero changes the optimal map. More details and references 
and given in section 1.5 of [BF012]. 

The discrete approximation of the combined problem (MA, BV2, C) in 
the generality of (1) is our present topic. 



4 JEAN-DAVID BENAMOU, BRITTANY D. FROESE, AND ADAM M. OBERMAN 

1.2. Contributions of this work. In order to enforce the state constraint 
(BV2) for all possible target sets Y, we use the signed distance function, 
d, to the boundary of the domain Y. Replace (BV2) by a Hamilton-Jacobi 
equation on the boundary: 

(HJ) H(Vu(x)) = d(Vu(x),Y) = 0, for x G dX. 

From the convexity of the domain, Y, the signed distance function, d, is also 
convex. This implies convexity of H. Using convexity of the solution, u, we 
show that (MA, BV2, C) is equivalent to (MA, HJ, C). However, the latter 
system is more amenable to computations, because it involves a local PDE 
(a convex Hamilton-Jacobi equation) instead of a global condition on the 
mapping. In full generality, any strictly convex Hamiltonian function which 
vanishes on dY could be used. This corresponds to the notion of defining 
function introduced in [Del91, Urb97]. The choice of the signed distance 
function is made for simplicity. 

We recall it has become common practice to use a distance function to 
determine a set, as is the case in the level set method. In that case, a 
Hamilton-Jacobi equation is used to solve for the distance function. What 
we are doing here is the converse, using the distance function to a particular 
set to represent the nonlinear PDE operator. 

Using this reformulation, we are then able to build approximations and 
prove convergence of the approximations in the setting of viscosity solutions 
(which include the regular case). Our main result, Theorem 3.6 below, is 
restated here. 

Theorem (Convergence). Let u be the unique convex viscosity solution 
of (MA, BV2, C). Suppose that py is Lipschitz continuous with convex sup- 
port, as in the assumption on the densities (1). Let u h,dS,da be a solution of 
the finite difference scheme (19). Then u h > ,da converges uniformly to u as 
h, d6, da — > 0. 

Remark. This theorem does not give a rate of convergence, which is typical 
of this kind of convergence result, since viscosity solutions can be singular. 
However, the filtered scheme which we use is formally second order accurate, 
so, in the case of smooth solutions, we expect second order accuracy. 

1.3. Weak solutions of the OT problem. Any convergent approach to 
the solution of the OT problem in the continuous setting must use an ap- 
propriate notion of solution. The optimal transportation problem with qua- 
dratic cost has a solution under quite general conditions on the measures to 
be transported. It frequently leads to singular solutions of (MA, BV2, C), so 
some notion of weak solutions is needed in the general setting. It is natural 
to ask under what conditions the system (MA, BV2, C) is well-posed. 

For the convenience of the reader, we briefly discuss regularity and weak 
solutions, for the purpose of using these notions of solutions for building 
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approximations. This is a substantial topic, which we can not cover here, 
we refer to [Vil03, Chapter 4] and the book [GutOl]. 

The existence of classical solutions to (MA, BV2, C) requires a smooth, 
convex target, and smooth strictly positive densities. The regularity theory 
for classical solutions can be found in [Del91, Urb97] and [Caf96]. 

Aleksandrov solutions [GutOl] are defined in terms of the Hessian measure 
of the potential function. This allows for the target measure to be singular 
with respect to Lebesgue measure. A typical example is the potential func- 
tion of a cone u{x) = ||x||, which is the potential for the mapping between 
a Dirac measure and a (scaled) characteristic function of the unit ball. 

Pogorelov solutions [Pog94] use specific target density in the form of a sum 
of weighted Dirac measures. The convex potential function is constructed 
by lifting affine functions with gradients prescribed by the support of the 
Dirac measures in the image. Adjusting functions changes the measure of 
the support of the affine function so that it matches the target value. 

Convex duality via the Legendre transform is the theoretical idea which 
links two notions of weak solution provided by Aleksandrov solutions and 
Pogorelov solutions. The basic building blocks are polyhedral functions, and 
their duals under the Legendre transform. 

Brenier's solutions [Bre91] are the most general discussed herein. They 
allow for the example of tearing a disc, see subsection 4.2. 

Viscosity solutions [Vil03, page 130], which require continuous density 
functions, are defined using the comparison principle. In the case where 
both measures are absolutely continuous, viscosity solutions coincide with 
Alexandrov solutions. 

For data in the form of Dirac masses, the constructive method of Pogorelov 
solutions is natural for the OT problem. In this setting, (BV2) is naturally 
satisfied, and the gradient mapping is approximated instead of the poten- 
tial, which reduced the discretization error of the mapping. This approach 
formed the basis of several early numerical methods: [OP88, CNP91, GM96] 
and more recently [BoslO, Merll]. Applying this approach to more general 
density functions requires quantization which introduces errors which are 
difficult to estimate. Also the best algorithm (assuming the initial density 
is also in the form of a sum of Dirac masses), the auction algorithm [Ber92], 
scales as 0(N 2 logN) in the number of Dirac masses. 

Viscosity solutions, though less general than Aleksandrov solutions, are 
well-suited to finite difference methods. They allow for optimal maps with 
discontinuous gradients to be computed. Densities are easily discretized on 
grids. Also, there is a robust convergence theory [BS91] which after the 
reformulation of (MA, BV2, C) into (MA, HJ, C) can be applied. Finally, 
in previous work we have developed fast and robust numerical methods for 
computing viscosity solutions of the Monge- Ampere equation with standard 
boundary conditions (see subsection 1.4). Moreover, the numerical results 
in [BF012] indicate that this approach reduces to the computational cost of 
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a linear elliptic solver which is log-linear with respect to the number of grid 
points. 

1.4. Relation to our previous work. This article builds on a series of 
papers which have developed solution methods for the Monge- Ampere equa- 
tion. The definition of elliptic difference schemes was presented in [Obe06], 
which laid the foundation for the schemes that followed. 

The first convergent scheme for the Monge-Ampere equation was built 
in [Obe08bj; it was restricted to two dimensions and to a slow iterative 
solver. Implicit solution methods were first developed in [BFO10], where it 
was demonstrated that the use of non-convergent schemes led to slow solvers 
for singular solutions. In [FOlla] a new discretization was presented, which 
generalized to three and higher dimensions; this also led to a fast Newton 
solver. 

The wide stencil schemes used in the convergent discretizations intro- 
duced a new parameter, the directional resolution, which led to decreased 
accuracy. In an attempt to improve accuracy on less singular solutions, a 
hybrid discretization was built in [FOllb]. This discretization combined 
the advantages of accuracy in smooth regions, and robustness (convergence 
and stability) near singularities. However, this was accomplished at the ex- 
pense of a convergence proof. In addition, it required a priori knowledge of 
singularities, which is not available in the OT setting. 

In [F012], the convergence theory of Barles and Souganidis was extended 
to filtered (nearly monotone) schemes. The filtered schemes replaced the 
hybrid schemes as a method to obtain more accuracy, without losing the 
convergence proof. On less singular solutions, the filtered schemes can be 
used with a compact stencil eliminated the complication and slightly in- 
creased cost of the wide stencil. 

The numerical resolution of (MA, BV2, C) was first addressed in [Frol2]. 
The method consisted of iteratively solving (MA) with Neumann boundary 
conditions, and projecting the resulting set onto the target set Y. The 
new projection is then used to derive new Neumann boundary conditions. 
This method required several iterations to reach a solution satisfying the 
state constraint (BV2). No convergence proof was available. We show in 
[BF012] that this method can be interpreted as a particular solution method 
for (MA, BV2, C). 

2. Representation and approximation of the boundary 

conditions 

In this section we describe our approach to approximating the boundary 
condition (HJ) based on the signed distance function to the boundary of the 
target set Y. 
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2.1. Representation and properties of the distance function. For a 

bounded, convex set Y, the signed distance function, 



(4) d(y,Y) 



is convex. 
The function 



+ dist(y, dY), y outside of Y, 
— dist(y, dY), y inside Y. 



H(y) = d(y,Y), 

which is now interpreted as a convex Hamiltonian, can be written in terms 
of the supporting hyperplanes to the convex set, 

(5) H(y) = sup {n(y ) ■ (y - y )} , 

yoedY 

where n(yo) is the outward normal to dY at yo; see Figure 1. Equivalently, 
since the image of the normals to dY is the unit sphere, we can write 

(6) H(y) = sup {n- (y-y(n))}, 

\\n\\=l 

where y(n) is the point in dY with normal n. The representations (5) and (6) 
are equivalent via duality: in the first case the convex set is represented by 
normals, in the second by points. We found the second representation to be 
useful, since sets are naturally approximated by a characteristic functions: 
this was the representation used in [BF012]. 

The statements above follow from the Supporting Hyperplane Theorem [BV04, 
Section 2.5], which says that if y$ £ dY, for a convex set Y, then yo has a 
(possibly non-unique) supporting hyperplane, 

P = {A(y) = 0\A(y) = n-(y-y )}, 

where A(y) < 0, for y & Y. Without loss of generality, ||n|| = 1, and we can 
define n to be (an) outward normal to Y at yo. Then, we can define 

H*(n) = n ■ y(n) = sup n • y, 
y&dY 

where the equality follows from the supporting hyperplane result. We have 
proven the following lemma. 

Lemma 2.1. Let H(y) = d(y,Y) be the signed distance function to the 
smooth, convex, bounded set, Y . Then 



H(y) = sup {y ■ n — H*(n)} 
IMI=l 



where 



(7) 



H*(n) = sup {y ■ n}. 

yoGdY 
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Figure 1. Polyhedral target set. 

Remark (Approximation of H*). Later we will make an approximation by 
taking the supremum over a finite subset of the admissible directions. These 
direction vectors are typically given by a uniform discretization of [0, 2tt], 
with angle discretization parameter da. We require only that da — > for 
convergence. 

2.2. Obliqueness. We recall here a fundamental property of maps char- 
acterized as the gradient of a convex potential. In [Del91, Urb97], this 
obliqueness result is used to prove existence of classical solutions to (MA- 
BV2). 

This condition, which leads to Lemma 2.3, will allow us to build an explicit 
monotone upwind discretization of (5), using points inside the domain. 




Figure 2. Illustration of the mapping y = Vtt(x) and the 
normal vectors. 
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Lemma 2.2. Suppose X is a convex domain, and Y = Vrt(X) is the im- 
age of X under the mapping Vu, where u is a convex twice continuously 
differ entiable function. Then the normal vectors make an acute angle, 

(8) n x -n y > 0. 
See Figure 2. 

Proof. Let Vu{x) = y G 8Y and let H(y) = dist(y, dY) be the signed 
distance function to Y, so that 

Y = {y | H(y) = 0}, 

and 

X = {x\ H(Vu(x)) = 0}. 
Then n y = VH(y) and, by the chain rule for differentiation, 

n x = cVD 2 u(x)VH(Vu(x)) = cD 2 u(x)VH(y) 

for a normalization constant c. Thus 

n x -n y = c{VH{y)) T D 2 u(x)VH(y) > 0, 

since convexity of u means D 2 u is positive definite. □ 

We are now able to give the characterization of (HJ) which will be used 
for the monotone discretization. 

Lemma 2.3. Let u G C 2 (X) be any convex function that satisfies (BV2). 
For any x £ dX with unit outward normal n x , the supremum in (HJ) can 
be restricted to vectors making an acute angle with n x : 

(9) H(Vu(x)) = sup {Vu(x) ■ n - H*(n) \ n ■ n x > 0} = 0. 

Hl=i 

Proof. The supremum above will be attained for a value of n = n* , which 
will be identical to the unit outward normal to the target at the point Vu(x). 

From Lemma 2.2, we know that n* ■ n x = n y ■ n x > 0. Consequently, it 
is only necessary to check values of n that make an acute angle with the 
boundary of the domain. □ 

2.3. Monotone discretization of H. In this section we explain how to 
build a monotone discretization of H using points at the boundary and on 
the inside of the domain X. 

The expression (9), which comes from writing the convex set Y in terms 
of its tangent hyperplanes, leads to a natural convergent finite difference 
discretization. This expression can be interpreted as a Hamilton-Jacobi- 
Bellman equation arising from an optimal control problem. 

We recall (see Definition 3.3) that an elliptic (monotone) discretization 
of H(Vu) takes a particular form. In the case at hand, the discretization 
is given in the following form. At the point X{ the discretization is a non- 
decreasing function of the differences between u(xi) and u(xj) where Xj are 
the neighbours of xi. 
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If we let T Xi = {n \ n ■ n Xi > 0, ||n|| = 1}, a simple way of writing an 
upwind discretization is to approximate the signed distance function by 

H(Vu( Xi )) = sup {Vu{ Xi )-n-H*{n)} « sup \ < x i)-< x i~ nh ) _ 
ner Xi n&i I h 

where h is a small discretization parameter. The discretization of the linear 
advection above can be directly implemented in the wide stencil framework 
as in [Obe08b] (see also section 3.3 below). As long as the domain is uni- 
formly convex and h is sufficiently small, obliqueness ensures that the point 
Xi — nh is inside the domain. This guarantees that the monotonicity for all 
n G Tj. Taking the supremum of monotone terms results in a monotone 
expression. 

Alternately, we can use standard compact upwind finite differences on 
a grid for the linear advection equation. Assume for simplicity that the 
boundary is linear and the local coordinate system such that at the boundary 
point xij, the normal is n Xij = (—1,0). Along this side, the set of admissible 
directions T Xi ■ will be given by 

{n = (ni,ri2) | n\ < 0, ||n|| = 1}. 

Then, letting dx denote the spatial resolution of the grid, we can approxi- 
mate the advection terms in (9) by 

v^.). BRiBl !M^) 

ax 



+ max{ri2 , 0} — h min{n2 , 0} - 



dx dx 

Due to the obliqueness property again, upwinding leads to a scheme that 
relies on values inside the square and, because n\ < 0, it is monotone. Taking 
the supremum of these monotone schemes over all admissible directions, we 
again preserve monotonicity of the scheme. 

We follow this last approach in [BF012] where we always use an initial 
density, px , whose support is embedded in a rectangular domain and padded 
by zeros. The robustness of our Monge- Ampere solver to degenerate solution 
allows for this simplification. 



3. Convergence 

We begin with a review of background material that will be needed to 
construct and prove the convergence of our scheme for solving the second 
boundary value problem for the Monge- Ampere equation. 

The viscosity approximation theory developed by Barles and Sougani- 
dis [BS91] provides criteria for the convergence of approximation schemes: 
schemes that are consistent, monotone, and stable converge to the unique 
viscosity solution of a degenerate elliptic equation. This framework was 
extended in [F012] to nearly monotone schemes, which are more accurate. 
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(11) F(x,u,p,M) 



3.1. Writing the combined problem as a single operator. It is con- 
venient to write the combined PDE and boundary conditions (MA, BV2, C) 
as a single (discontinuous) operator on X 

(10) F(x, u(x),Vu(x),V 2 u(x)) = 0, x G X 

where F depends on px,PY and H: 

f det(M) - px(x)/py(p), x£X 
\H(p), x G dx, 

along with the convexity condition (C). 

Example 3.1. The Dirichlet problem is simply given by replacing H(p) with 
u — g for x G dX. 

3.2. Convergence theory. In this subsection we review the definitions and 
the general theory for convergence schemes, it is applied to the equation (11) 
in a subsequent section. 

Definition 3.2 (Consistent). The scheme F e is consistent with the equa- 
tion F = if for any smooth function <f> and x G f2, 

limsup F*(y, <f>(y) + £, </>(•) + £) < F*(x, 4>(x),V(j)(x), D 2 (j)(x)), 
e— >0,y— >x,£— >0 

liminf F\y, <j>(y) + £, 0(-) + > F*(x, «/»(x), V^(x), ,D 2 0(x)). 
e-+0,j/— >:r,§— >0 

Definition 3.3 (Elliptic). T/ie scheme F € is elliptic if it can be written 

F € [v] = F e (x, v{x),v(x) - v(-)), 
where F e is nondecreasing in its second and third arguments, i. e. 

(12) s<t, «(•)<«(■) F e (x,s,u(-)) < F £ (x,t,u(-)) 

Definition 3.4 (Nearly Elliptic). The scheme F € is nearly elliptic if it can 

be written as 

(13) F e [v] = F M [u] + Fp[u) 

where Fm is a monotone (elliptic) scheme and Fp is a perturbation, which 
satisfies 

lim llFpll = 0. 

Using these definitions, we now recall the main convergence theorem 
from [F012], which is an extension of the convergence theory of [BS91]. 
This almost monotone scheme is related to the work of Abgrall [Abg09] 
which is a "blending" of a monotone and a higher-order scheme. 

Theorem 3.5 (Convergence of Nearly Monotone Schemes [F012]). Let u be 

the unique viscosity solution of the PDE (10) and let u e be a stable solution 
of the consistent, nearly elliptic approximation scheme (13). Then 

u e — > u, locally uniformly, as e — > 0. 



12 JEAN-DAVID BENAMOU, BRITTANY D. FROESE, AND ADAM M. OBERMAN 



Moreover, if the non-monotone perturbation Fp is continuous, u e exists and 
is stable. 

The filtered approximation schemes combine a monotone (elliptic) scheme Fm 
with a more accurate, non- monotone scheme Fa- These schemes, which are 
a particular type of nearly monotone scheme (thus they are convergent by 
the theorem above), have the form 

'F A [u]-Fl,[u\ 



(14) 



F*[u] = F M [u]+r(e)S 



re 



where r(e) — > as e — > 0. Here the function S, which is called a filter, can 
be defined, for example, by 



(15) 

See Figure 3. 



S(x) 



x \\x\\ < 1 

[|ac|| > 2 

-x + 2 l<x<2 

-x-2 -2<x<-l. 




Figure 3. Filter function 



3.3. Discretization of Monge- Ampere equation. The equation we want 
to solve is 

det(D 2 u(x)) = p x {x)/p Y {Vu[x)) + (u). 
We first describe the elliptic (monotone) scheme for the Monge-Ampere 
operator, which underlies the filtered scheme. This scheme was developed 
in [FOlla, Frol2]. 

Remark. It is somewhat complicated to describe the discretization. First, we 
combine a discretization of the Monge-Ampere operator with the convexity 
constraint. In addition, regularization terms are added to make the operator 
differentiable when Newton's method is applied. Further, the dependence 
on the gradient in (MA) means that small changes in the values of u can 
lead to large changes in the equation. This causes instabilities when noncon- 
vergent schemes are used. We also modify the discretized Monge-Ampere 



NUMERICAL MONGE-AMPERE FOR OPTIMAL TRANSPORTATION 13 



operator to compensate for the dependence on the gradient. We describe 
the modifications step by step. 

To begin, we use Hadamard's inequality to represent the determinant of a 
positive definite matrix, det(M) < Himu, with equality when M is diagonal. 
Then, we can write 

det(M) = min{ni(0 T MO)ii | T = 1} 

where O is an orthogonal matrix. This last inequality, applied to the Hessian 
of a convex function, corresponds to taking products of second derivatives 
of the function along orthogonal directions, 

d 

det(D 2 u)= min J] u Ujl/j , 
{ui...u d }ev 

where V is the set of all orthonormal bases for 

In the special case where the source density px vanishes, the Monge- 
Ampere operator reduces to the convex envelope operator [Obe07, OS09] 
which corresponds to directional convexity in each direction. The convex 
envelope operator is approximated by enforcing directional convexity in grid 
directions [Obe08b, Obe08a]. 

In the current setting, the convexity constraint is enforced by additionally 
replacing the directional derivatives with their positive part. In addition, 
to prevent non-convex solutions when the right-hand side vanishes, we will 
also subtract the negative parts of these second derivatives, 

{d d 
IT max{u u . u . , 0} + mm{u„ v , 0} 

which is valid when u is convex. These modifications ensure that a non- 
convex function cannot solve our Monge- Ampere equation (with non-negative 
right-hand side). 

We discretize the operator above in two ways. First, we make the conven- 
tional step of replacing derivatives by finite differences. Second, instead of 
considering all orientations, we replace V with a finite subset which uses only 
a finite number of vectors v that lie on the grid and have a fixed maximum 
length. The second discretization is called the the directional discretization, 
and we quantify it using d9, the angular resolution d6 of our stencil (see 
Figure 4). In this figure, values on the boundary are used to maintain the 
directional resolution d9 (at the expense of lower order accuracy in space be- 
cause the distances from the reference point are not equal). Another option 
is to narrower stencils as the boundary is approached, which leads to lower 
angular resolution, but better spatial resolution. We denote the resulting 
set of orthogonal vectors by Q. 
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(a) (b) (c) 

Figure 4. Neighboring grid points used for width one 
(green), two (yellow), and three (blue) stencils. The illus- 
tration shows the neighbors in the first quadrant. The mod- 
ification near the boundary is illustrated in the second and 
third figures. 



Each of the directional derivatives in the Monge- Ampere operator is dis- 
cretized using centered differences: 

T> vv Ui = n * „ (u(xi + vh) + u(xi - vh) - 2u(xi)) . 
\\ u \\ h 

In order to handle non-constant densities, we also need to discretize the 
gradient. This is accomplished using the Lax-Friedrichs scheme, which uses 
centered differences augmented by a small multiple of the Laplacian to en- 
sure monotonicity. 

d 

-px/pviyu) « -px/py (T> Xl u, . . .,V Xd u) + 5^V, 

0=1 



XjXjU- 



The centered difference discretization of the first derivatives is given by 
V v Ui = — (u(xi + vh) - u(xi - vh)) . 

To preserve monotonicity, we require the parameter 5 to satisfy 

(16) 5 > Kh, 

where K is the Lipschitz constant (with respect to y) of px{x) / py{v)- 
Then an elliptic discretization of the Monge- Ampere equation is 

(17) MA h f> S [u]= min Gf* [,] 



NUMERICAL MONGE-AMPERE FOR OPTIMAL TRANSPORTATION 15 



where each of the G^' de ' S , \u] is defined as 
(ui,...,v d ) l 

d d 

G tn' 5 ^d) M = II max{D„ J „ ;I .«, 0} + ^ mii^XV^u, 0} 

(18) 

- px(x)/p Y (J) Xl u, . . .,V Xd u) + 8y)T> XiXi u - u(x ). 

This monotone scheme forms the basis of the filtered scheme (14). For 
improved accuracy on smooth solutions, we combine the monotone scheme 
with the accurate scheme Fa, which is simply a standard centered difference 
discretization of the (two-dimensional) equation 

u XlXl u X2X2 -u 2 XlX2 - p x (x)/p Y (u Xl ,u X2 ) -u(x ). 

We denote the resulting discretization by MAs[u]. A similar discretization 
is easily constructed in higher dimensions by using standard finite differ- 
ence discretizations for the other terms. Additional details can be found 
in [FOllb]. 

3.4. Convergence to the viscosity solution. We combine the almost 
monotone schemes for (MA, C) with the upwind, monotone scheme for (HJ) 
into one equation, which we show converges to the unique convex viscosity 
solution of the system (MA, BV2, C). The combined scheme is given as 

(19) F h,de,da [u] . _ \F h / e W\ : ,,eX 



J I 



H h > da \u\i xi G dX. 



In this definition, Fp is the filtered scheme for the Monge-Ampere equa- 
tion (14), which relies on the discretizations described in subsection 3.3, 
and H is the upwind discretization of the boundary condition described in 
subsection 2.3. 

Theorem 3.6 (Convergence). Let u be the unique convex viscosity solution 
of (MA, BV2, C). Suppose that py is Lipschitz continuous with convex sup- 
port, as in the assumption on the densities (1). Let u h) ,da be a solution of 
the finite difference scheme (19). Then u h > dd)da converges uniformly to u as 
h, d9, da —> 0. 

Proof. From Theorem 3.5, we need only verify that the scheme is consistent 
and nearly elliptic (nearly monotone). 

Consistency and near monotonicity of the scheme for the Monge-Ampere 
equation have been established in [FOlla, Frol2]. 

We recall the form of the boundary condition in (9), 

H(Vu(x)) = sup {Vu(x) ■ n — H*(n) \ n • n x > 0}. 



This is discretized using forward or backward differences for the gradient, 
which are consistent as h — > 0. The supremum is further approximated 
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by restricting to a finite subset of directions, with angular resolution da. 
Since the Legendre-Fenchel transform H*(n) is continuous, and since these 
admissible directions are approximated with an accuracy on the order of da, 
this approximation is consistent as da — > 0. 

By exploiting the obliqueness property (Lemma 2.3), we were able to con- 
struct an upwind discretization of the boundary condition, which is mono- 
tone by construction. □ 



We now provide some brief computational results to demonstrate that 
the approximation scheme described in this paper can be used in practice 
to numerically solve the optimal transportation problem. Full details of 
the numerical method and its implementation are given in a companion 
paper [BF012]. 

For the first two examples, which admit exact analytical solutions, we 
provide the maximum norm of the distance between the mappings obtained 
from the exact and computed solutions (u ex and u comp respectively): 



In the tables, the number of grid points used to approximate each di- 
mension of the domain X and target Y are proportional to Nx and Ny 
respectively. 

Remark (Reading mappings from the figures (5-6) ). The source density is 
embeded in a squared domain discretized by a cartesian grid and padded 
with zeros. We represent on the left the restriction of the grid to the sup- 
port of the source and on the right , its image by the optimal map. The 
optimal mapping can be interpreted from the figures by noting, first, that 
the centre of masses are mapped to each other. Next, moving along a grid 
line in the source, the corresponding point in the target can be found by 
using monotonicity of the map: the corresponding grid point is in the same 
direction. 

4.1. Mapping an ellipse to an ellipse. We consider the problem of map- 
ping an ellipse onto an ellipse. To describe the ellipses, we let M x ,M y be 
symmetric positive definite matrices and let B\ be the unit ball in M. d . Now 
we take X = M X B\, Y = M y B2 to be ellipses with constant densities /, g 
in each ellipse. 

In M 2 , the optimal map, which is linear, can be obtained explicitly [MO04] . 
It is given by 



4. Computational Results 



max ||Vti ea; (x) — Vit, 



•comp{.%) || 2- 



Vu(x) = M y R e M~ x 




cos(#) 
sin(0) 



-sin(0) 
cos(0) 
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the angle is given by 

tan(0) = trace(M 2 T 1 M~ x J) / trace^" 1 M~ l ) 
and the matrix J is equal to 



J = R 



it/2 



-1 

1 



We use the particular example 

0.8 
0.4 



M x 



0.6 0.2 
0.2 0.8 



which is pictured in Figure 5. 

Results are presented in Table 1, which demonstrates first order accuracy 
in both Nx and Ny- 



1 

0.5 


-0.5 




-0.5 0.5 1 

(a) 



1 

0.5 


-0.5 




Figure 5. (a) An ellipse X and (b) its image under the 
gradient map Vit (§4.1). 









Maximum Error 






Nx 






N Y 








8 


16 


32 


64 


128 


256 


32 


0.1163 


0.0773 


0.0693 


0.0669 


0.0665 


0.0062 


64 


0.1188 


0.0403 


0.0302 


0.0291 


0.0282 


0.0283 


128 


0.1214 


0.0302 


0.0201 


0.0174 


0.0168 


0.0168 


256 


0.1206 


0.0278 


0.0116 


0.0101 


0.0092 


0.0091 


362 


0.1175 


0.0291 


0.0098 


0.0063 


0.0057 


0.0056 



Table 1. Distance between exact and computed gradient 
maps for map from an ellipse to an ellipse. 
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4.2. Mapping from a disconnected region. We now consider the prob- 
lem of mapping the two half-circles 

X = {(31,2:2) I xi < -0.2, (xi + 0.2) 2 + x\ < 0.85 2 } 

U {(xi,x 2 ) I xi > 0.1, (xi - 0.1) 2 + x\ < 0.85 2 } 

onto the circle 

Y = {( yi ,y 2 ) I yl + yl < 0.85 2 }. 
This example, which is pictured in Figure 6, is degenerate since the domain 
is not simply connected. Nevertheless, our method correctly computes the 
optimal map, as the results of Table 2 verify. 









Maximum Error 






N x 






Ny 








8 


16 


32 


64 


128 


256 


32 


0.0453 


0.0267 


0.0255 


0.0258 


0.0259 


0.0258 


64 


0.0397 


0.0184 


0.0158 


0.0146 


0.0144 


0.0139 


128 


0.0392 


0.0097 


0.0063 


0.0066 


0.0065 


0.0064 


256 


0.0432 


0.0110 


0.0084 


0.0087 


0.0086 


0.0073 


362 


0.0448 


0.0130 


0.0070 


0.0047 


0.0045 


0.0039 



Table 2. Distance between exact and computed gradient 
maps for map from two semi-circles to a circle. 



4.3. Other geometries. To give a flavor of the types of geometry that 
can be captured by solving (HJ) with our approximation scheme, we also 
provide several different computed maps in Figure 7. These were all obtained 
by mapping a square with uniform density onto a specified convex set, whose 
boundary could consist of a combination of straight and curved edges. While 
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no exact solutions are available, we can certainly verify that the computed 
gradient does map the square into the correct target set. 

5. Conclusions 

A rigorous approach to finding solutions of the Optimal Transportation 
problem was presented, using convergent finite difference approximations to 
the elliptic Monge- Ampere PDE with corresponding boundary conditions. 

The target domain was represented using the signed distance function 
to the boundary. This representation led to a Hamilton-Jacobi equation 
for the optimal transportation boundary conditions. This Hamilton-Jacobi 
equation was discretized using a monotone finite difference scheme, which 
relied on obliqueness conditions for the solution to build a scheme using only 
values of the function inside the source domain. 

The combined PDE and boundary conditions led to an nearly elliptic 
finite difference scheme, which was proved to converge to the unique viscosity 
solution of the underlying PDE. The convergence proof used the theory of 
filtered schemes, [F012], which is an extension of the monotone schemes in 
the Barles-Souganidis theory [BS91]. 

Computational results were presented which are consistent with the pre- 
dictions of the theory. A complete numerical implementation of the scheme, 
including a fast Newton solver, is presented in the companion paper [BF012]. 

A first application is to use this method to compute JKO gradient flows [JK098, 
OttOl]. Another application is to OT problem with more general cost func- 
tions c{x — y), using a version of the Monge- Ampere PDE for c-convex 
potentials. 
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